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Abstract 

LQG  compensator  control  of  transverse  vibrations  was  implemented  on  an  aluminum  can¬ 
tilevered  beam  in  a  ’’smart  structure”  paradigm.  The  beam  was  mounted  with  two  self-sensing, 
self-actuating  piezoceramic  patches.  The  Euler- Bernoulli  beam  equation  was  discretized  via  a 
Galerkin  type  approximation  (referred  to  as  the  full  order  model).  To  reduce  the  size  of  the  re¬ 
sulting  Unite  dimensional  approximating  system,  the  Proper  Orthogonal  Decomposition  (POD) 
was  employed  as  a  reduced  basis  method.  A  reduction  of  dimension  from  34  to  2  was  obtained 
through  the  model  reduction  technique.  Feedback  control  based  on  the  reduced  order  system  was 
implemented  in  real  time  using  a  dSpace  DS1103  control  system.  Experimental  results  indicate 
that  POD  based  control  achieves  comparable  control  attenuation  with  full  order  model  based 
control. 


1  Introduction 

Real-time  control  of  smart  material  structures  requires  numerical  integration  of  the  linear 
system  at  each  discrete  time  step  at  which  the  real-time  processor  runs.  For  complex 
structural  systems  such  as  shells  [16,  17,  18]  or  structural  acoustic  systems  [6],  the  compu¬ 
tationally  intensive  algorithms  needed  for  the  online  computation  of  PDE-based  controls 
must  be  supplemented  by  model  reduction  techniques,  Numerical  studies  in  [5]  have  shown 
the  feasibility  of  employing  model  reduction  techniques  in  feedback  control  of  thin  cylin¬ 
drical  shells  perturbed  by  a  periodic  external  force.  This  investigation  provides  the  next 
step  in  the  development  of  efficient  model  based  feedback  control  methodologies  by  im¬ 
plementing  the  ideas  presented  in  [5]  in  a  physical  system.  We  chose  to  demonstrate  the 
ideas  in  the  context  of  a  beam  system  but  the  extension  of  the  ideas  to  more  complex 
systems  is  rather  straightforward. 

The  proliferation  of  the  use  of  the  Proper  Orthogonal  Decomposition  (POD)  as  a 
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reduced  basis  method  in  control  and  estimation  applications  can  be  attributed  to  its 
ability  to  accurately  represent  the  system  data  with  only  a  small  number  of  basis  elements 
(or  POD  modes).  The  system  data  obtained  at  different  parameter  values  can  either 
be  experimental  measurements  of  the  state  or  numerical  simulations  of  the  system.  The 
POD  method  attempts  to  extract  characteristic  information  from  the  data  through  an 
orthogonal izal  ion  procedure.  It  is  essentially  a  linear  transformation  of  a  multivariate 
data  set  into  an  optimal  set  of  uncorrelated  variables  or  POD  modes.  Approximation 
applications  of  POD  include  turbulent  coherent  flows  [3,  12,  13,  27,  38,  20],  structures 
[30,  31,  32,  33,  35],  materials  processing  [24,  23,  29,  28,  41],  characterization  of  human 
faces  [39],  and  statistical  and  pattern  recognition  fields  [19],  to  name  a  few.  In  numerical 
control  applications,  the  use  of  POD  in  open  loop  control  of  fluids  described  by  Navier 
Stokes  type  PDE  can  be  found  in  [29,  28]  while  its  first  use  in  feedback  control  design  (in 
the  context  of  smart  material  structures)  is  reported  in  [5].  Other  recent  applications  in 
feedback  control  can  be  found  in  [1,  2,  24,  23,  25]  while  POD  use  in  inverse  problems  for 
damage  detection  is  first  found  in  [9]. 

Another  type  of  reduced  basis  method,  the  Lagrange  model  reduction  technique,  has 
been  numerically  studied  in  [21,  22]  (in  the  context  of  approximation  and  control  of  fluid 
flows)  and  in  [5,  15]  (in  approximation  and  control  of  thin  shell  dynamics).  This  model 
reduction  technique  employs  the  snapshots  themselves  as  basis  functions.  As  detailed  in 
[5,  15],  the  method  possesses  inherent  difficulties  since  the  user  must  rely  on  knowledge  of 
the  solution  of  the  system  in  order  to  obtain  linearly  independent  snapshots.  The  question 
of  which  parameter  values  (or  time  instances  in  the  case  vibrating  structures)  to  use  for 
the  snapshots  in  order  to  capture  the  key  characteristics  of  the  system  is  another  difficult 
issue  to  be  addressed  in  the  Lagrange  reduced  basis  method.  These  questions,  as  we  shall 
see  in  Section  3.2,  are  of  little  concern  in  the  POD  approach. 

In  this  paper,  we  employ  the  POD  model  reduction  technique  in  real-time  experimental 
control  of  a  beam  structure.  The  simplicity  of  the  beam  structure  (i.e.,  the  relatively 
low  dimension  of  the  resulting  system  in  the  Galerkin  expansion  using  full  order  basis 
functions  such  as  cubic  splines)  allows  us  to  implement  real-time  full  order  based  control 
and  compare  these  results  with  control  based  on  the  POD  reduced  order  method.  In 
Section  2,  the  governing  equations  are  presented  and  existing  theories  regarding  well- 
posedness  of  the  model  are  summarized.  Approximation  methods  regarding  the  full  order 
method  and  POD  reduced  order  method  are  presented  in  Section  3.  Control  methodologies 
are  given  in  Section  4  and  experimental  results  can  be  found  in  Section  5.  Finally,  we  offer 
some  concluding  remarks  in  Section  6  relating  our  POD  based  methods  to  standard  finite 
element  and  modal  based  control  techniques. 


2  Beam  Model 

In  this  section,  we  briefly  present  the  equations  governing  transverse  displacements  of  an 
Euler- Bernoulli  beam  with  Kelvin- Voigt  damping  and  cantilevered  boundary  conditions. 
We  also  summarize  existing  results  regarding  well-posedness  of  the  model.  We  denote  the 
beam  length,  width,  thickness,  density,  Young’s  modulus  and  Kelvin- Voigt  damping  by 
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Figure  1:  Cantilever  beam  with  piezoceramic  patches 


£,b,h,  p,E  and  cd,  respectively.  The  origin  is  taken  to  be  at  the  clamped  edge  of  the 
beam,  and  the  axial  direction  is  denoted  by  the  .x-axis  (see  Figure  1).  A  pair  of  identical 
piezoceramic  patches  are  bonded  on  opposite  sides  of  the  beam  with  edges  located  at  x-\  and 
.x 2 .  Passive  patch  contributions  arising  from  material  changes  due  to  the  presence  of  the 
patches  are  included  in  the  model.  Patch  parameters  are  denoted  by  the  subscript  pe,  thus, 
the  patch  thickness,  width,  density,  Young’s  modulus  and  Kelvin- Voigt  damping  are  given 
by  hpe,  bpe ,  ppe,  Epe  and  C£>pe}  respectively.  We  denote  the  piezoelectric  constant  relating 
mechanical  strain  and  applied  electric  held  by  d3 j.  Finally,  we  denote  the  transverse 
displacement  by  y  and  the  voltages  applied  to  the  front  and  back  patches  by  V\  and  V2, 
respectively. 

As  derived  in  [11],  the  transverse  (or  bending)  equation  of  the  beam  is  given  in  terms 
of  resultant  moments  b Mx  by 

d2y  d2(bMx)  d2(bMx)pe 

Pdt 2  dx 2  +  dx2  1  '  U 

Passive  patch  contributions  are  incorporated  in  the  model  above  and  hence  the  linear  mass 
density  p(x)  =  phb  +  2 bppehpexpe{x)  is  piecewise  constant.  The  characteristic  function 
Xpe{x)  employed  to  isolate  patch  contributions  is  defined  by 


XPe{x) 


1  ,  x1<x<x2 
0  ,  otherwise  . 


(2) 


Incorporating  both  internal  damping  and  material  changes  due  to  the  presence  of  the 
patches,  the  internal  moment  resultant  bMx  has  the  form 


0*2 

bMx(t,x )  =  -£Y(x)—  y(t,x)  -  cDI(x ) 


d3 


dx2dt 


y(t , 


X 


(3) 
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where 


- —  h3b  2b 

EI(x)  =  E—  +  —Epea3xPe{x) 


—  h36  2b 

cDI(x)  =  cD—  +  — <i),. .  o-AXr'  |  (4) 


and  g3  =  ( h/2  +  hpe )3  —  ft3/ 8  (we  refer  the  reader  to  [11,  Chapter  3]  for  details  regarding 
patch  contributions  to  the  internal  moment  resultant).  When  voltages  are  applied  to  the 
front  and  back  patches,  the  induced  external  moment  ( bMx)pe  is  given  by 


(bMx)pe(t,x)  =  ^Epebd3i(h  +  hpe)xPe(a’)[Vi(f)  -  V2{t)] 


(5) 


External  transverse  forces  acting  on  the  beam  are  modeled  by  the  function  f(t,x).  Can¬ 
tilever  boundary  conditions  are  given  by 


=  0)=0, 

and  initial  conditions  are  denoted  by 

2/(0,  x)  =  y0{x) 


.!/.,(/,  I)  =  ^ f)  =  0  , 


|f(0,  *3  =  ■ 


(6) 


(7) 


To  reduce  regularity  requirements  in  the  discontinuous  moment  resultants  we  employ  the 
weak  formulation  of  (1)  given  by 

/"  ;  w  ,  /•'  -  ,  d2y(t,x) 

J0  Phb — Qp — 4>(x)dx  +  J  2bf)pehpe — — — <p{x)dx 

,  ->>  ,,  d2y(t,x)d2(f)(x) 

dx  T  /  . ,  Eped 3 

«/07l  O 


+ 


l 


+ 


cU2  dx 2 

1  d3y(t,  x)  d24>(x) 


I  cDI- 
J  0 


dx2dt  dx 2 
’^2  1 


dx-\- 

J  X 


X1 

2b 
x\  3 

2b 

Xl 


dx2  dx2 

d3y(t,x)d2<f>{x) 

jcmPed3  dx2gt 

d2cj)(x ) 


dx 


dx 


(8) 


/  ^M3a(ft  +  hpe)(y1(f)  -  v2 (t))^-£p-dx 


+ 


fi 

/  f{t,x}<j)(x)dx  , 
Jo 


for  all  </(x)  G  H  =  ft/|(0,£)  =  {</>  G  ift2(0,  £)|</>(0)  =  <//(0)  =  0}.  Here  /  =  6ft3/ 12,  the  state 
space  is  taken  to  be  the  Hilbert  space  H  =  i2(0,£),  and  we  seek  solutions  y(t ,  ■)  G  H. 

Using  a  variational  framework,  the  existence  and  uniqueness  of  the  solution  to  (8) 
has  been  established  in  [11,  Chapter  4]  and  those  results  will  be  summarized  here.  We 
reformulate  (8)  in  a  variational  setting  by  first  defining  the  sesquilinear  forms  a \  and  a2 
defined  on  V  X  V  -*•  d'  by 


<7 1 


(4Vf|i=  and  cr2((/>,  ip)  =  (cDI(f)",  4>") 


(9) 
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for  all  (p,if)  G  V.  Note  that  V  and  H  form  a  Gelfand  triple  V  H  ~  H*  <— >  V* , 
i.e. ,  U  is  continuously  and  densely  embedded  in  H  and  we  identify  H*  with  H  through 
the  Riesz  map  (see  [42,  pages  165,261]).  We  take  the  duality  pairing  (vfy*  y  to  be  the 
extension  by  continuity  of  the  inner  product  (•,  •) H  from  V  X  H  to  V*  X  H .  We  then  define 
iht +  f  eV*  by  (( bMx)pe  +  f)  ((f))  =  (( bMx)pe(t ),  </>")  +  </(/),  </>)  for  <f>  G  V.  Equation 
(8)  can  now  be  abstractly  formulated  as 

#yqy  +  °2(y(t),  4>)  +  4>)  =  (( bMx)pe(t )  +  /(t),  </>)y%y  .  (10) 

It  can  be  readily  verified  that  (i)  the  stiffness  sesquilinear  form  0\  is  symmetric,  V- 
continuous  and  V -elliptic,  and  (ii)  the  damping  sesquilinear  form  a2  is  U-continuouS  and 
y-elliptic.  Under  the  additional  assumptions  (iii)  the  forcing  terms  (bMx)pe  +  /  satisfy 
the  regularity  condition  ( bMx)pe  +  /  £  U2((0,  T),  W),  (iv)  y0  £  V  and  (v)  y1  G  H  then 
Theorem  4.1  in  [11]  (which  is  first  proven  in  [8])  states  that  there  exists  a  unique  solu¬ 
tion  y  of  (10)  with  the  regularity  properties  y  (E  X2((0,  T),  V),  y  G  h2((0,T),U)  and  y  G 
T2((0,  T),  V*).  Moreover,  the  solutions  depend  continuously  on  the  data  (y0:  y i,  (bMx)pe  + 
/)  in  that  the  map  (y0,  y i,  (bMx)pe  +  /)  — >  ( y ,  y)  is  continuous  from  V  X  H  X  i2((0,  T),  \  “j 
to  T2((0,  T),V)  x  T2((0,  T),V). 

For  control  applications,  the  equivalent  operator  formulation  of  the  abstract  variational 
form  (10)  is  more  familiar.  Again  we  give  a  brief  summary  of  the  theoretical  results.  Under 
the  assumption  that  o\  and  <r2  are  U-continuous,  we  can  define  the  operators  Ai,  A2  from 
V  —)■  V*  given  by 


)(<p)  ••  a,(c).  rj  v  (A2if))(<f>)  =  (r2(<f>,if))  ,  (11) 


for  all  (f),  if>  G  V.  Then  an  equivalent  formulation  of  equation  (10)  is  given  by 

py(t)  +  A\n(l.)  +  A2y(l.)  =  ( bMx)pe(t )  +  f(t)  in  V*  .  (12) 

The  hrst  order  formulation  of  (12)  is  obtained  by  first  defining  the  product  spaces  V  = 
V  x  V  and  7i  —  V  X  H .  We  then  let 


rin  = 


'  0 

„  1  , 

o 

o 

.  fit)  . 

5  *3  —  —Epebdsi{h  -)-  hpe) 

.  Xpe{x)  Xpe(x )  . 

(13) 


and  u(t)  =  [  V i  ( /. ) .  V2(t)]T .  Finally  we  dehne  the  system  operator  A  in  matrix  form  by 

^4  = 


0  I 
— A\  —  A2 

doinM  =  [x  =  ((f),  (p)  G  'H\ip  G  V  and  A\<p  +  ^2^  G  H} 
and  the  strong  hrst  order  form  of  (12)  is  now  given  by 


(14) 


w(t)  =  Aw(t)  +  Bu(t)  +  A(t) 

w(0)  =  w0  =  (y0,yi)  , 


(15) 
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where  w  =  (y,y).  The  mild  solution  to  (15)  is  given  by 


w(t)  —  T(t)w o+  f  T(t  —  s)[Bu(s)  +  lF(s)\ds  (16) 

Jo 

and  its  existence,  uniqueness  and  equivalence  to  the  variational  solution  is  detailed  in 
[11].  Here  T  is  the  semigroup  generated  by  A  (see  [34]  for  details  regarding  semigroup 
generators). 


3  Galerkin  Discretization  and  Model  Reduction 

A  Galerkin  method  employing  cubic  splines  as  basis  functions  was  used  to  discretize  the 
system  (8).  As  discussed  in  [16],  cubic  splines  were  chosen  over  modes  since  for  the  sys¬ 
tem  under  consideration,  explicit  expressions  of  the  modes  and  natural  frequencies  are  not 
readily  available.  In  order  to  use  mode  shapes  in  the  expansion,  they  must  first  be  numeri¬ 
cally  approximated  or  experimentally  determined.  Moreover,  the  piecewise  discontinuities 
in  the  density,  stiffness  and  damping  coefficients  arising  due  to  contributions  by  the  piezo- 
ceramic  patches  cause  difficulties  for  the  use  of  standard  finite  element  packages  since 
the  meshes  must  be  aligned  with  the  regions  bounded  by  the  patches.  Moreover,  finite- 
difference  approximations  are  not  employed  due  to  difficulties  associated  with  high-order 
equal  ions  in  the  beam  model. 

The  standard  Galerkin  approximation  employing  cubic  splines  is  efficient  for  numeri¬ 
cal  simulations  but  the  size  of  the  resulting  finite  dimensional  system  is  still  not  suitable 
for  real-time  feedback  control  implementation.  Hence  we  investigate  the  use  of  Proper 
Orthogonal  Decomposition  reduced  basis  methods  in  decreasing  the  dimension  of  the  ap¬ 
proximating  system. 

In  Section  3.1,  we  briefly  discuss  the  Galerkin  discretization  method,  and  specify  the 
cubic  spline  basis  functions  employed.  We  denote  the  use  of  modified  cubic  spline  basis 
functions  in  the  discretization  by  full  order  methods.  In  Section  3.2,  the  ideas  underlying 
the  Proper  Orthogonal  Decomposition  reduced  basis  method  are  presented.  This  method 
is  referred  to  as  the  reduced  basis  method.  Finally,  in  Section  3.3,  we  report  in  detail 
the  methods  and  discretization  sizes  employed  in  numerically  implementing  the  methods 
presented  in  this  section. 

3.1  Galerkin  Approximation 

The  approximate  solution  y  £  L2((0,T),V)  to  (8)  is  obtained  by  replacing  the  space 
V  with  appropriate  finite  dimensional  subsets  =  span{HJ}^1  C  II  =  L2( 0,£)  and 
performing  the  Galerkin  expansion 

yM(t,  x)  =  ■  (17) 

4  =  1 

The  approximating  finite  dimensional  subsets  V  'V  are  chosen  to  satisfy  the  condition 
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V 


0  as 


(HI)  for  any  (f>  £  V,  there  exists  a  sequence  ^  in  V  A  such  that  <f>^  —  (f> 

Af  — *  oo . 

We  then  restrict  the  weak  form  (8)  to  V  'A'  with  basis  functions  also  used  as  test  functions 
and  hxed.  (All  of  the  resulting  matrices  are  understood  to  depend  on  J\f  throughout, 
but  for  convenience  we  shall  suppress  this  in  our  notation.)  This  yields  the  matrix  system 
approximating  (8)  given  by 


{M  +  Mpe )  z{t)  +  (D  +  Dpe )  z(t )  +  (  K  E  +  K  Epe)  z(t )  —  A(^)  T  Bu(t ) 
2(0)  =  z0,  z( 0)  =  z%  , 


(18) 


where  z(t)  =  [z^t), .  .  . ,  zj\f[t)]T  is  the  vector  of  coefficients.  The  matrices  in  (18)  are 
defined  by 


[M]k,i  =  phb  Bi(x)Bk{x)dx 

[KE]kj  =  El  f  B'/(x)B''(x)d 
Jo 

[D]k,i  —  cDI  (  B"  (x)B"(x)dx 
Jo 

[B]k(t)  =  f  f(x,t)Bk(x)dx 
Jo 

u(t)  =  [V4(t) ,  v2{t)Y 


f-X  2 

[Mpe]  k,l  ^bppehpe  /  Bi(x)Bk(x)dx 

J  X\ 

[ KEpe ]kJ  =  ^ Epea3  r  B\!{r)B'l\x)d 

O  J  x\ 

2  b  rx? 

[DPe]k,i  =  — cDpea3  /  B"(x)B"(x)dx 

O  J  x  1 

[B]k,  1  =  lEpebd31(h  +  hpe)  f  B”(x)d . 

J  X-l 


[&\k,2  —  —\Epebd31{h  +  h 


pe 


^2 

xi 


Bk(x)d 


x  . 


The  first  order  reformulation  of  (18), 

w(t)  =  Aw{t )  +  Bu(t)  +  F(t) 
w( 0)  =  w0  =  [z0,  Zjf  , 

is  obtained  by  letting  w  —  (z,z)T, 


m  = 


A  = 

■  I 

0 

0  M  +  Mj 

r  _ 

„  "1 

-1  ■ 

I 

0 

|_  0  M  +  Mpe 

-1 


0 

F(t) 


0  I 

—  {KE  +  KEpe)  -( D  +  Dpe ) 


,  and  B  = 


■  I 

0 

-1 

'  0 

0 

M  +  Mpe 

B 

(19) 


(20) 


(21) 

(22) 


Cubic  splines  were  employed  as  basis  functions  in  the  standard  (or  full  order)  Galerkin 
expansion.  This  choice  of  basis  functions  was  motivated  by  smoothness  requirements,  ac¬ 
curacy,  adaptability  to  different  boundary  conditions  and  flexibility  with  regard  to  internal 
and  external  patch  contributions.  To  define  the  standard  cubic  splines,  we  first  partition 
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the  interval  [0,  l]  with  grid  points  xn  =  nhx ,  hx  =  £/N,  n  =  0, .  .  .  ,  N .  The  standard  cubic 
splines  are  then  defined  by 


f  -  a’n-2)3, 


bn(x)  = 


hi 


G  [Xyy  —  2  7  Xn_ l] 


hi  +  3 h2x(x  -  Xn- 1)  +  3ha-(x  -  J'n-l)2  -  3(2;  -  Xn-2)3,  X  G  [xn-l,  Xn] 
h3x  +  3/i2(xK+i  -  x)  +  3ha.(2:n+i  -  x)2  -  3(s„+i  -  x)3,  x  G  [xn, xn+i] 


(xn+2  -  x)3, 

0, 


X  G  [Xn-|-1  ,  Xyy-^] 
otherwise  , 


(23) 

for  n  =  —  1, .  .  . ,  N  +  l,(see  [36]).  These  standard  splines  must  be  modified  to  satisfy  the 
clamped  edge  boundary  conditions  (6).  The  modified  cubic  splines  satisfying  the  zero  slope 
and  zero  displacements  conditions  at  x  —  0  are  given  by 

f  h0(x)  —  2h_i(x )  —  26i(x)  ,  n  =  1 

BJx)  =  {  ;  n  ;  y  J  ’  24 

■  l  bn{x)  ,  n  =  2,  •  ■  • ,  JV  +  1  ,  k  J 

for  a  total  of  N  +  1  basis  functions.  These  basis  elements  satisfy 


Bn(0)  =  B'n(0)  =  0  n  =  l,---,N  +  l. 

Therefore,  if  N  standard  cubic  splines  are  employed,  the  system  (20)  is  of  dimension 
M  —  2 (TV  +  1)  (note  that  the  factor  2  comes  from  reformulating  the  system  in  first  order 
form). 

A  general  framework  for  the  convergence  of  full  order  solutions  has  been  given  in 
the  electronic  version  of  [7].  Under  the  assumption  that  (HI)  is  satisfied  by  the  cubic 
splines,  Lemma  4.1  in  [7]  guarantees  convergence  of  the  full  order  solution  to  the  infinite 
dimensional  solution  of  (8). 


3.2  Basis  Reduction  via  Proper  Orthogonal  Decomposition 

The  POD  reduced  order  basis  method  we  discuss  here  involves  the  use  of  a  smaller  approx¬ 
imating  subspace  of  .  To  differentiate  full  order  basis  functions  from  reduced  basis 

elements,  we  denote  the  ith  full  order  basis  functions  by  Bi  while  we  denote  the  ith  POD 
reduced  basis  function  by  $j.  We  also  denote  the  full  order  dimension  index  by  M  and 
reduced  order  dimension  index  by  Mp. 

To  create  the  POD  basis  elements,  we  take  Ns  temporal  snapshots  of  the  model, 
,  x) }  ,  and  seek  basis  elements  of  the  form 

N‘ 

=  J2  of v* (*}-•> x)  ■  (25) 

3= 1 

The  snapshots  can  either  be  numerical  solutions  or  physical  measurements  of  the  state.  We 
require  each  basis  element  to  resemble  all  the  snapshots  in  the  sense  that  it  maximizes 

1  m ..... 

s  3  =  1 


,  subject  to  =  ||$j||2  =  1  .  (26) 
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As  detailed  in  [5,  12,  15,  29],  the  coefficients  of  the  ith  POD  basis  function  (he.,  the  a*-  in 
(25))  are  the  elements  of  the  ith  eigenvector  of  the  covariance  matrix  C,  where  C  is  defined 
by  . 

[C]M=  Y  (yN\tk, .<■),  !!*(<<,•))  ,  k,e=l,...,N,  .  (27) 

Since  C  is  Hermitian  and  nonnegative,  we  can  order  the  eigenvalues  Ai  >  A2  >  .  .  .  > 
Ajvs  A  0  and  it  has  a  complete  set  of  orthogonal  vectors.  Note  that  the  eigenvectors  are 
normalized  in  such  a  way  that 


h  / 

Ol  •  Oi  —  < 


0,  k  f  £ 

1 


l  NsXk  ’ 


k  =  £ 


The  POD  basis  functions  <I>,  were  shown  to  be  orthonormal  in  [29].  It  was  also  shown  in 
[12]  that  the  numerical  solution  using  the  POD  basis  functions  in  the  expansion  is  optimal 
in  the  sense  that  it  maximizes  the  time  average  of  the  displacements.  This  optimality  of 
the  POD  basis  functions  is  summarized  in  [12,  Prop.  2.3]. 


Lemma  1  Let  {$i,  $2,  •  •  • ,  $jv«}  denote  the  orthonormal  set  of  POD  basis  elements  and 
Ai  >  .  .  .  >  Ajva  denote  the  corresponding  set  of  eigenvalues.  If  yNs  =  ^=4  bi(t)<$>i  denotes 
the  approximation  to  y  with  respect  to  this  basis,  then  for  any  arbitrary  orthonormal  basis 
{ipi,ip2,  ■  ■  ■ ,  ipN,} ,  the  following  hold 


1.  <  bi(t)b*{t )  =  A i&ij,  where  <C  ■  denotes  the  average  over  time, 

Ns  Ns  Ns 

2.  for  every  Ns,  ^  <C  bi(t)b*(t )  >>  =  Ea8>E  <  at(t)a*(t)  >  ; 

where  the  ai(t )  are  coefficients  of  the  approximation  to  y ,  y ^  u,(/.)( using  the 

arbitrary  orthonormal  basis  }  E  • 

Lemma  1  provides  a  systematic  way  of  creating  an  Mp  dimensional  POD  reduced 
basis  (where  jVp  <C  Ns )  out  of  Ns  snapshots.  Since  the  most  significant  POD  modes 
correspond  to  the  largest  eigenvalues,  we  take  jVp  to  be  the  smallest  integer  such  that 
Yfpfi  A i  —  E(=1  Aj.  The  ratio  gives  the  percentage  of  the  mean  square'  error  when  the  first 
Mp  POD  modes  are  eliminated  in  the  reduced  basis  representation  (see  [24]  for  details  and 
the  references  cited  therein).  Initial  approximation  results  for  the  use  of  POD  elements  in 
Galerkin  schemes  are  given  in  [26]. 

As  emphasized  in  [5],  we  note  that  for  the  POD  reduced  basis  method,  a  fundamental 
practical  question  is  not  of  convergence  but  of  (i)  how  to  take  an  ensemble  of  snapshots 
which  lead  to  a  good  approximation  of  the  full  order  solution  and  (ii)  how  to  choose  a 
reduced  set  of  basis  elements  that  will  capture  the  approximation  properties  of  the  original 
ensemble.  It  is  difficult  to  address  the  first  question  since  it  requires  knowledge  of  the 
solution  of  the  original  physical  system.  Since  the  POD  reduced  basis  method  does  not 
use  the  snapshots  as  basis  themselves,  the  question  of  independence  of  snapshots  is  not 
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important.  Hence  one  way  to  capture  the  characteristics  of  the  full  order  solution  is  to 
take  a  sufficiently  large  number  of  snapshots.  This  is  an  acceptable  strategy  since  the 
covariant  matrix  C  and  the  POD  basis  functions  are  computed  offline.  In  relation  to  the 
second  question,  Lemma  1  provides  a  systematic  method  to  choose  a  reduced  set  of  basis 
elements  out  of  the  Ns  snapshots.  The  lemma  shows  how  well  the  first  Mp  POD  basis 
functions  approximate  the  displacement  average  contained  in  all  the  snapshots  (see  [5]). 


3.3  Numerical  Implementation 

Numerical  discretization  using  the  full  order  method  was  performed  using  16  standard 
cubic  splines,  resulting  in  an  M  =  34  dimensional  matrix  system  (20).  Integration  of  the 
basis  functions  (see  (19))  were  performed  using  a  four  point  Gaussian  quadrature  rule  with 
Nq  =  64  nodes.  This  quadrature  is  of  the  form 


N„  4 


[  f(x)dx  «  ckf{xkl )  , 


(28) 


1=1  k= 1 


where  the  quadrature  weights  c*  and  points  x^i  on  the  Ith  interval  [/  —  1,  /]  *  h„ ,  hq  =  £/Nq 
are  given  by 


Cl  = 


c2  = 


c3  = 


c4  = 


49 

hq 

6(18  +  V^0) 

2 

49 

hq 

6(18  -  a/30) 

2 

49 

hq 

6(18  -  V30) 

'  y 

49 

hq 

6(18  +  V^0) 

2 

Xll  =  hq 

X21  =  hq 
X3 1  =  hq 

X41  —  hq 


1  ^5  +  2^30 

2  2^35 

1  y/l5  -  2^30 

2  2y^5 


1  y  15  —  2\/30 

2+  2^35 


1  \/ 15  +  2  \/30 

2+  2^35 


(see  [40]  for  details). 

To  obtain  the  snapshots  for  creating  the  POD  basis  elements,  the  uncontrolled  full 
order  system  excited  by  a  voltage  spike  to  the  back  patch  was  numerically  simulated  on 
the  interval  t  =  0  to  t  =  0.5  seconds.  The  voltage  spike  is  triangular  in  shape  with 
a  duration  of  0.001  .s  and  a  magnitude  at  the  peak  of  90  V.  A  variable  stepsize  stiff 
ODE  solver  was  employed  in  time  stepping  the  system.  One  hundred  evenly  spaced  time 
snapshots  ( Ay  =  100)  were  taken  from  which  the  100  X  100  covariant  matrix  C  given 
by  (27)  was  created.  After  obtaining  the  eigenvectors  and  eigenvalues  of  the  covariant 
matrix,  the  number  of  POD  basis  functions  was  determined  by  ensuring  the  percentage 
of  the  displacement  average  of  all  snapshots  contained  by  the  reduced  basis  was  at  least 
90%.  It  can  be  seen  in  Table  1  that  this  requirement  is  satisfied  even  with  1  POD  basis 
function. 
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Afp 

i 

2 

3 

4 

5 

6 

.9008399 

.9999942 

.9999988 

.9999996 

1.0 

1.0 

Table  1:  Ratio  of  full  order  model  snapshots  captured  by  the  POD  basis  functions;  Ns  = 
100  snapshots;  Afp  is  the  number  of  POD  reduced  basis  elements. 


4  LQG  Compensator  Control 

Infinite  dimensional  control  for  the  system  (15)  where  T  =  0  or  T  is  periodic  has  been 
studied  extensively  in  [7,  11,  37].  Since  our  goal  is  to  implement  real  time  control,  we 
will  concentrate  on  finite  dimensional  control  methods  for  the  matrix  equation  (20)  and 
consider  only  transient  shell  vibrations  (F  =  0).  This  case  we  consider  is  suitable  when 
controlling  vibrations  caused  by  an  impulse  force  (such  as  a  hammer  impact),  or  when  the 
structure  responds  to  a  perturbation  and  vibrates  to  a  steady  state. 

A  state  estimator  must  be  employed  in  the  control  scheme  due  to  the  nature  of  the 
sensing  device  employed  in  the  experiment.  Vibrations  measured  at  only  one  point  x  on 
the  beam  by  the  proximity  probe  are  used  as  observations.  Thus,  we  consider  observations 
of  the  state  given  by  the  operator  C  :  V  — >  Y.  wob  =  Cw(t),  where  Y  is  the  observation 
space  and  wob(t )  =  y^(t:x)  ~  //(/.,  ;r). 

We  present  the  approximate  LQG  control  problem  with  no  exogenous  input  in  Sec¬ 
tion  4.1.  In  order  to  experimentally  implement  the  control  method  in  Section  4.1,  we 
must  modify  the  assumptions  that  observations  can  be  obtained  at  any  time  through¬ 
out  the  time  interval  of  interest,  and  that  the  state  and  state  estimator  can  be  exactly 
integrated.  Thus,  we  discuss  the  discrete-time  compensator  problem  in  Section  4.2. 


4.1  Continuous-time  Finite  Dimensional  Control 

The  infinite  horizon  control  problem  is 

roo 

minimize  .J(u,w0)  =  /  {(Qw(t),w(t))  (Ru(t),u(t))}dt 

Jo 

subject  to  the  matrix  equation 

w(t)  =  Aw(t )  +  Bu{t )  ,  u>(0)  =  w0 

wob(t)  --  Cw(t )  . 


(29) 


(30) 


Here  C  is  the  observation  matrix  which  maps  w( /)  with  the  observation  wob(t ),  while  Q 
and  R  are  control  parameter  matrices  used  to  weigh  the  state  and  control,  respectively.  We 
follow  the  control  methods  numerically  implemented  in  [5,  18]  where  the  optimal  control 
is  given  by 

u{t)  =  — Kwc(t )  .  (31) 
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Here  the  state  estimator  is  given  by 


wc(t)  =  Awc(t)  +  Bu(t)  +  F[w0b(t)  —  Cwc(t)] 
Wc( 0)  =  wCo  =  Cw0  , 


(32) 


(see  [11]  for  details)  and  the  optimal  feedback  gain  (regulator)  K  and  observer  (Kalman 
filter)  gain  F  are  given  by 

K  =  R  '/h'  ll 

(33) 

f  =  n^R-1 , 

where  n  and  n  are  solutions  to  the  algebraic  Riccati  equations 

n  a  +  at  n  -  nBR-'B^n  +  q  =  o 
flAT  +  AII-  \  \CtR  '<+jfl  +  0  =  0, 


(34) 


respectively.  Playing  a  similar  role  to  that  of  Q  and  R  in  the  control  problem,  the  matrices 
Q  and  R  are  design  criteria  for  the  state  estimator. 

The  structure  of  the  observation  matrix  C  depends  on  the  sensor  employed  in  the 
experiment.  In  our  model,  the  sensor  is  a  proximity  probe  located  at  the  back  of  the  beam 
sensing  displacements  at  the  point  x.  The  observation  matrix  C  is  thus  of  the  form 


C 


Bi(xr%..*.,Bjj-(x),  0,  ■  ■  ■  ,0  , 

AT  . 


where  the  /£,(+)+  are  the  basis  functions  evaluated  at  x. 

It  was  shown  in  [7,  11]  (  see  also  [4])  that  the  approximate  controls  and  Riccati  so¬ 
lutions  based  on  full  order  methods  converge  (as  M  — >  oo)  to  the  corresponding  infinite 
dimensional  optimal  control  and  Riccati  operators.  The  assumptions  needed  to  establish 
convergence  of  the  full  order  Riccati  operators  and  optimal  controls  are  that  (i)  (HI)  holds, 
(ii)  V  is  compactly  embedded  in  H ,  and  (iii)  the  infinite  and  finite  dimensional  systems  are 
stabilizable  and  detectable.  For  the  POD  based  controls,  the  questions  of  convergence  are 
not  meaningful.  The  issue,  rather,  is  how  the  reduced  order  based  controls  would  perform 
when  implemented  in  physical  infinite  dimensional  systems.  In  a  previous  investigation, 
the  issue  was  computationally  addressed  when  POD  reduced  order  based  controls  were 
applied  to  the  full  order  system  in  controlling  thin  cylindrical  shell  vibrations  (see  [5]).  In 
Section  5,  we  further  address  the  issue  of  effectiveness  on  the  orginal  infinite  dimensional 
system  by  implementing  POD  based  controls  on  a  physical  system. 


4.2  Discrete-Time  Finite  Dimensional  Control 

In  experimental  applications,  the  signals  from  the  sensors  are  digitized  and  the  real-time 
processor  can  only  perform  at  a  discrete  sample  rate  At.  Thus,  the  state  estimator  equation 
(32)  can  only  be  evolved  in  time  in  discrete  time  steps,  and  the  control  voltage  can  be 
only  computed  at  this  rate.  The  numerical  ODE  approximation  method  to  solve  the 
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state  estimator  equation  must  satisfy  the  following  criteria:  (i)  the  control  u(tj)  must  be 
calculated  before  the  arrival  of  the  data  at  the  next  time  step  /,J+1  =  tj  +  At,  (ii)  the 
method  must  be  sufficiently  accurate  to  resolve  system  dynamics  and  (iii)  since  the  ode 
systems  are  often  stiff,  the  method  must  be  .4-stable  or  a-stable.  We  chose  a  modified 
backward  Euler  method  given  in  [11,  Chapter  8.2.1],  The  fast  sample  rate  (and  hence 
small  At)  at  which  we  can  carry  out  the  experiment  allows  the  use  of  this  method.  An 
.4-stable  modified  backward  Euler  method  integrating  the  state  estimator  (32)  at  time  t j+i 
is  given  by 

wCj+1  =  (I  -  AtAc)-^wCj +  At(I  -  AtAjj-'Fwobftj) 

(35) 

=  R(Ac)wCj  +  A tR{Ac)Fwoh{tj)  , 

where  Ac  =  A  —  BK  —  FC ,  R[AC)  =  (/  —  At.4c)_1  and  the  constant  time  step  is  At  = 
tj+ 1  —  tj.  Note  that  the  method  is  modified  from  standard  backward  Euler  methods  since 
the  observation  wob{tj)  at  future  time  steps  are  not  available.  We  now  summarize  the 
discrete-time  algorithm  in  Algorithm  1,  which  is  essentially  Algorithm  8.5  in  [11] 


Offline 

0) 

(ii) 

(iii) 

(iv) 

Construct  matrices  .4,  B ,  (7,  Q ,  R ,  Q ,  R 

Solve  Riccati  equations  (34)  for  11  and  11 
Construct  K  =  R~x BTU,  F  =  \  \CT R 
and  Ac  =  A  —  BK  -  FC 

Construct  R(AC)  =  (I  —  AtAc)-1  and 
R{AC)F  =  (I-  AtAJ-'F 

Online 

(i) 

Collect  observation  w0b(tj ) 

(ii) 

Time  step  the  discrete  compensator  system 

=  R{Ac)wCj  +  AtR(Ac)Fwob(tj) 

(iii) 

Calculate  the  voltage  u(tj+ 1)  =  —KwCj+1 

Algorithm  1:  Discrete  compensator  for  systems  running  at  a  sample  rate  of  At. 


5  Real-Time  Feedback  Control  Implementation 

In  Table  2,  we  report  the  dimensions  and  parameters  of  our  experimental  beam  struc¬ 
ture  depicted  in  Figure  2.  The  aluminum  beam  parameters  p ,  E,  c d  and  the  lead  zirconate 
titanate  piezoceramic  patch  parameters  ppe,  Epe,  ci)v, .  d31  were  obtained  from  the  manufac¬ 
turers.  Typically,  parameter  estimation  techniques  must  be  carried  out  in  order  to  obtain 
more  accurate  values  but  since  the  goal  of  this  investigation  is  the  application  of  the  POD 
reduced  basis  method  in  real  time  control,  we  do  not  consider  the  inverse  problem  in  detail 
here.  Rather,  we  refer  the  reader  to  [11,  Chapter  5.4]  for  discussions  and  experimental 
results  regarding  estimation  of  parameters  and  inverse  problems  involving  the  beam. 
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Beam 

Patch 

£  =  0.286  m 
h  =  0.001  m 
b  =  0.2543  m 
p  =  3.438  x  103  kg/m 3 

E  =  7.062  x  1010  N/m 2 
cD  =  1.04  X  106  Ns/m 2 
x  =  0.11076  m 

hpe  —  5.3  x  10-4  m 
ppe  =  7.45  x  103  kg/m 3 

Epe  =  6.4  x  1010  N/m 2 
cDpe  —  3-96  x  10s  Ns/m 2 
d3 1  =  262  x  10-12  m/V 
x-i  =  0.02041  m 
x-2  =  0.04592  m 

Table  2:  Beam  and  Patch  Parameters 


Figure  2:  Experimental  beam  with  piezoceramic  patches 


Numerical  simulations  were  performed  to  obtain  reasonable  values  of  the  control  pa¬ 
rameter  matrices  Q.  li.  Q  and  R  to  be  used  in  (29)  and  (34).  We  sought  parameters  leading 
to  maximum  control  voltages  within  the  ilOOE  range  of  the  patches  while  at  the  same 
time  providing  good  attenuation  (see  [5,  17]  for  discussions  regarding  the  choice  of  control 
parameters  in  shell  control  simulations).  The  matrices  employed  were  of  the  form 


Q  —  d\ 

Q  =  di 


KE  +  KEpe  0 

0  M  +  Mpe 

IN  0 
0  IN 


R  =  rjm  ,p  =  1 
R  =  rjsxs  ,5  =  1, 


(36) 


where  p  is  the  number  of  actuators  and  s  is  the  number  of  sensors.  In  Table  3,  we  report 
the  values  of  d1:  rj,  eh,  for  three  discretization  sizes:  full  order  system,  POD  system 
with  2  POD  basis  functions,  and  POD  system  with  1  POD  basis  function.  To  illustrate 
the  expected  performance  of  the  control,  we  also  report  in  Table  3  the  maximum  of  the 


14 


real  parts  of  the  eigenvalues  of  the  matrices  A,  A  —  BK  and  A  —  FC  for  each  of  the  three 
discretization  sizes.  It  can  be  seen  that  the  maximum  real  part  of  the  eigenvalues  of  the 
closed  loop  system  A—BK  was  3  orders  of  magnitude  less  than  that  of  the  system  matrix  A, 
while  for  the  state  estimator,  A  —  FC  was  2  orders  of  magnitudes  less.  To  further  illustrate 
the  features  of  the  reduced  basis  method,  the  controllability  matrices  of  the  linear  systems 
arising  from  the  three  discretization  sizes  were  analyzed.  Controllability  is  a  necessary  (but 
not  sufficient)  condition  for  a  linear  system  to  be  driven  to  zero  starting  at  any  initial  state. 
Linear  systems  described  by  the  matrices  A  and  B  in  the  first  order  linear  system  (20)  are 
controllable  if  the  2Afpx2sAfp  controllability  matrix  C(A,  B)  =  [_B|A_B|A2i?|  ■  •  ■  \A2J^P~1B], 
where  A  =  A2J^P ,  B  =  B2J^P ,  has  full  rank  2 Afp.  For  systems  with  rank  (C(A,  B ))  =  m.  < 
2 jVp,  any  initial  condition  of  the  state  z2^p  in  the  m  dimensional  subspace  of  II { 2J^P 
spanned  by  the  columns  of  C(A,B)  can  be  driven  to  zero  (see  [14]).  Generally  speaking, 
the  larger  the  deficit  of  C(A,B)  from  its  maximum  possible  rank,  the  more  difficulty  in 
stabilizing  or  controlling  the  linear  system  (30).  The  ranks  of  the  controllability  matrices 
for  the  three  discretization  systems  are  reported  in  the  last  row  of  Table  3.  It  can  be  seen 
that  the  full  order  linear  system  is  more  difficult  to  stabilize  than  that  of  the  two  POD 
systems  (whose  controllability  matrices  are  of  full  rank). 


Full  Order 

£ 

II 

to 

Afp  =  1 

<h ,  n 

di,h 

max(3i(eig(A))) 

2  x  108,  0.98 

1  x  103, 1 
— 3.107e  -  2 

6  x  1010, 1  x  1(T4 

5  x  10s,  1  x  KT1 
— 3.111e  -  2 

IBB— 

max(3i(eig(A  —  BK))) 

—  1.215e  +  1 

— 6.625e  +  1 

— 6.695e  +  1 

max(:R(eig(/l  —  FC))) 

— 8.526e  +  0 

— 3.445e  +  0 

— 3.433e  +  0 

rank(6'(  A,  B)) 

9/34 

4/4 

2/2 

max  possible  rank 

Table  3:  Control  Parameter  Values  and  some  properties  of  the  resulting  linear  systems 
due  to  the  three  discretization  sizes:  full  order,  POD  system  with  2  basis  functions  and 
POD  system  with  1  basis  function. 

In  Figure  3,  we  present  a  diagram  of  the  experimental  setup  and  implementation  of 
the  online  component  of  Algorithm  1.  Voltage  spikes  to  the  back  patch  (to  excite  the 
beam)  were  generated  by  a  DS1103  dSpace  control  system.  The  excitation  signal  was  low 
pass  filtered  and  amplified  before  being  applied  to  the  back  patch.  The  voltage  spike  was 
amplified  so  as  to  produce  90  volts  at  the  peak.  A  proximity  probe  located  at  the  back  of 
the  beam  at  x  —  x  —  0.11076  m  was  used  to  measure  displacements  and  the  observation 
readings  were  digitized  through  one  analog  to  digital  channel  of  the  dSpace  hardware.  This 
observation  signal  enters  the  online  component  of  Algorithm  1  as  w0b(tj).  By  employing  the 
discrete  modified  backward  Euler  method  (35),  the  state  estimator  wc{t1+ 1)  was  obtained 
and  multipled  with  the  gain  matrix  K  to  produce  the  control  voltage.  The  control  signal 
was  then  low  pass  filtered  and  amplified  before  being  sent  to  the  front  patch.  A  constant 
discrete  time  rate  of  A t  =  10-4  s  was  employed  in  running  the  real-time  processor. 
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dSpace  ADC  dSpace  DAC  dSpace  DAC 

Channel  17  Channel  1  Channel  2 


Figure  3:  Experimental  setup  and  implementation  of  online  component  of  Algorithm  1. 

In  Figure  4,  we  report  the  uncontrolled  and  controlled  displacements  in  the  left  column 
and  the  control  voltages  in  the  right  column  for  the  three  discretization  sizes:  full  order 
method,  POD  based  control  with  2  POD  basis  elements  and  POD  based  control  with  1 
POD  basis  element.  Note  that  the  three  control  systems  have  basically  attenuated  the 
displacements  after  1  second. 
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6  Concluding  Remarks 

Our  initial  experimental  results  suggest  that  application  of  POD  based  control  in  real 
time  control  is  most  promising.  Creation  of  the  POD  reduced  basis  using  numerical 
solutions  as  snapshots  is  relatively  easy  to  implement  since  a  lot  of  effort  has  been  invested 
in  developing  accurate  models  and  efficient  numerical  approximation  methods  for  smart 
material  structures.  The  online  component  control  method  (shown  in  Algorithm  1)  can  be 
efficiently  implemented  in  real-time  due  to  the  simplicity  of  the  algorithm  and  the  simple 
programming  of  the  real-time  control  implementation  system  employed. 

A  comparison  of  the  attenuation  levels  due  to  three  discretization  sizes,  i.e. ,  full  order 
(34  dimensional  linear  system),  POD  with  2  basis  functions  (4  dimesions)  and  POD  with 
1  basis  function  (2  dimensions)  illustrates  the  feasibility  of  POD  based  control  in  phys¬ 
ical  implementation.  Aside  from  decreasing  the  dimension  of  the  linear  approximating 
system,  the  POD  model  reduction  also  results  in  linear  systems  that  are  controllable  and 
observable. 

For  the  particular  physical  system  (a  cantilevered  beam)  used  here  to  illustrate  the 
effectiveness  of  POD  based  control  methods,  one  might  expect  that  an  effective  low  order 
model  (dimension  <  3)  would  suffice  based  on  the  research  literature  on  modal  based 
methods  in  structures.  However,  this  does  not  lessen  the  importance  of  our  contributions 
here.  First,  we  note  that  the  POD  based  methods  offer  an  attractive  implementable  al¬ 
ternative  to  modal  methods  since  analytical  modes  for  even  the  simple  piezo  mounted 
beam  are  not  attainable.  Moreover,  such  systems  do  not  have  normal  modes  [10],  further 
complicating  any  modal  based  computational  control  method.  Furthermore,  the  compu¬ 
tation  of  POD  basis  elements  is  even  simpler  than  computation  of  modal  elements.  When 
compared  to  standard  finite  element  approaches,  the  POD  based  methods  offer  control 
authority  with  much  smaller  approximating  systems  and  hence  are  clearly  preferable  for 
implementation.  In  addition,  for  more  complicated  structural  problems  such  as  those 
arising  in  many  applications,  the  POD  based  methods  will  offer  a  far  superior  approach 
in  implementation  of  control  and  stabilization  design.  Finally,  for  control  design  in  fluid 
or  electromagnetic  systems  where  modes  are  difficult  if  not  impossible  to  find  and  where 
finite  element  approximations  result  in  really  large  (W  on  the  order  of  10,000-100,000) 
approximating  systems,  POD  based  methods  offer  even  more  exciting  possibilities 
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Figure  4:  Displacements  at  x  =  0.11075  m\  (a)  Full  Order  (b)  2  POD  (c)  3  POD 
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